Ejemplo 2: Modelación de casos de COVID

Motivación

Al comenzar la pandemia producida por COVID-19, surgió el interés de analizar la distribución espacial de los riesgos infección.

Ajustamos un modelo lineal generalizado con efectos aleatorios espaciales con INLA para estimar riesgo de aparición del virus en diferentes radios censales de la Ciudad de Córdoba, Argentina.

Datos

Se trabajará con una base de datos vectorial (.gpkg) que incluye polígonos de los radios censales de Córdoba y Gran Córdoba. Para cada radio censal se cuenta con información de

  • Cantidad de habitantes en miles de personas (Poblacion)
  • Nivel de fragmentación urbano (fragmentacion)
  • Valor de la Tierra promedio
  • Porcentaje de hogares con nececsidades básicas insatisfechas (HogaresNBI)
  • Porcentaje de hogares con hacinamiento (HogariesHacinamiento)
  • Porcentaje de hogares con Jefe/Jefa con al menos nivel de estudio universitario (HogariesJefe.Univ.)
  • Porcentaje de hogares con computadoras (HogaresConComput.)
  • Cantidad de habitantes por \(km^2\) (HabSup)
  • Cantidad de bancos dentro del radio censal (Bancos)
  • Cantidad de farmacias dentro del radio censal (Farmacias)
  • Cantidad de centros de salud dentro del radio censal (EstSaludInt)
  • Cantidad de bancos, farmacias y centros de salud del radio censal (Ban_Farm_Salud)
  • Cantidad de casos detectados de COVID-19 (Casos)

Lectura base de datos

Carga paquetes
library(sf)
library(spdep)
library(tmap)
library(INLA)
datos <- st_read("data/Base_07_07_radios.gpkg", quiet = TRUE)
  • Cálculo de número de casos esperados por radio censal
datos$E <- datos$Poblacion * sum(datos$Casos) / sum(datos$Poblacion)
tm_shape(datos) +
  tm_polygons(fill = 'E')

  • Cálculo de cociente de infección estandarizada (Standardized Infection Ratio, SIR) por radio censal
datos$SIR <- datos$Casos / datos$E
tm_shape(datos) +
  tm_polygons(fill = 'SIR')

  • Referencia para el efecto aleatorio y el término del error
datos$re_u <- 1:nrow(datos)
datos$re_v <- 1:nrow(datos)

Ajuste del Modelo

Generación de lista con vecindarios

# Generación de lista con vecindarios
nb <- poly2nb(datos)

head(nb)

plot(st_geometry(datos), border = "grey", lwd = 0.5)
plot(
  nb,
  coordinates(as(datos, "Spatial")),
  add = TRUE,
  col = "blue",
  points = FALSE,
  lwd = 0.5
)

## [[1]]
##  [1]    2    3    4   29   37  237  251  252  258  266  267 1498 1500
## 
## [[2]]
## [1]  1  3  8 14 37
## 
## [[3]]
## [1] 1 2 4 8
## 
## [[4]]
## [1]   1   3   5   6   7   8 266 267
## 
## [[5]]
## [1]   4   6 267
## 
## [[6]]
## [1]   4   5   7  10  11 270
# Generación de vecindarios para INLA
file_adj <- tempfile("map.adj1")
nb2INLA(file_adj, nb)
g <- inla.read.graph(filename = file_adj)
# Ajuste del Modelo inflado en ceros
summary(datos)


formula  <-
  Casos ~ 
  Fragmentacion +  
  HogaresNBI + 
  HogaresHacinamiento +
  HogaresJefe.Univ. + 
  Bancos + 
  f(re_u, model = "besag", graph = g) + 
  f(re_v, model = "iid")


res_inflpoi <-
  inla(
    formula,
    family = "zeroinflatedpoisson1",
    data = datos,
    E = E,
    control.predictor = list(compute = TRUE)
  )
summary(res_inflpoi)
##    Poblacion      Fragmentacion    ValorTierra       HogaresNBI     
##  Min.   :   7.0   Min.   :1.000   Min.   :   200   Min.   : 0.0000  
##  1st Qu.: 629.8   1st Qu.:2.000   1st Qu.:  2450   1st Qu.: 0.3755  
##  Median : 835.0   Median :4.000   Median :  6000   Median : 1.0855  
##  Mean   : 873.2   Mean   :3.367   Mean   : 17074   Mean   : 1.6939  
##  3rd Qu.:1086.0   3rd Qu.:4.000   3rd Qu.: 12500   3rd Qu.: 2.4040  
##  Max.   :3273.0   Max.   :4.000   Max.   :220000   Max.   :20.6897  
##  HogaresHacinamiento HogaresJefe.Univ. HogaresConComput.     HabSup       
##  Min.   :0.0000      Min.   : 0.0000   Min.   : 0.00     Min.   :  0.006  
##  1st Qu.:0.1364      1st Qu.: 0.9802   1st Qu.:11.96     1st Qu.: 44.619  
##  Median :0.5051      Median : 3.3943   Median :18.05     Median : 73.769  
##  Mean   :0.8617      Mean   : 5.5653   Mean   :20.02     Mean   : 83.815  
##  3rd Qu.:1.2639      3rd Qu.: 9.0126   3rd Qu.:25.26     3rd Qu.: 97.168  
##  Max.   :6.6606      Max.   :26.6055   Max.   :57.17     Max.   :606.024  
##      Bancos          Farmacias       EstSaludInt     Ban_Farm_Salud  
##  Min.   :0.00000   Min.   :0.0000   Min.   :0.0000   Min.   :0.0000  
##  1st Qu.:0.00000   1st Qu.:0.0000   1st Qu.:0.0000   1st Qu.:0.0000  
##  Median :0.00000   Median :0.0000   Median :0.0000   Median :0.0000  
##  Mean   :0.05301   Mean   :0.2446   Mean   :0.1536   Mean   :0.4512  
##  3rd Qu.:0.00000   3rd Qu.:0.0000   3rd Qu.:0.0000   3rd Qu.:1.0000  
##  Max.   :8.00000   Max.   :6.0000   Max.   :6.0000   Max.   :9.0000  
##      Casos                    geom            E                 SIR        
##  Min.   : 0.0000   MULTIPOLYGON :1660   Min.   :0.002241   Min.   : 0.000  
##  1st Qu.: 0.0000   epsg:22174   :   0   1st Qu.:0.201595   1st Qu.: 0.000  
##  Median : 0.0000   +proj=tmer...:   0   Median :0.267299   Median : 0.000  
##  Mean   : 0.2795                        Mean   :0.279518   Mean   : 1.007  
##  3rd Qu.: 0.0000                        3rd Qu.:0.347649   3rd Qu.: 0.000  
##  Max.   :36.0000                        Max.   :1.047749   Max.   :89.967  
##       re_u             re_v       
##  Min.   :   1.0   Min.   :   1.0  
##  1st Qu.: 415.8   1st Qu.: 415.8  
##  Median : 830.5   Median : 830.5  
##  Mean   : 830.5   Mean   : 830.5  
##  3rd Qu.:1245.2   3rd Qu.:1245.2  
##  Max.   :1660.0   Max.   :1660.0  
## 
## Call:
##    c("inla.core(formula = formula, family = family, contrasts = contrasts, 
##    ", " data = data, quantiles = quantiles, E = E, offset = offset, ", " 
##    scale = scale, weights = weights, Ntrials = Ntrials, strata = strata, 
##    ", " lp.scale = lp.scale, link.covariates = link.covariates, verbose = 
##    verbose, ", " lincomb = lincomb, selection = selection, control.compute 
##    = control.compute, ", " control.predictor = control.predictor, 
##    control.family = control.family, ", " control.inla = control.inla, 
##    control.fixed = control.fixed, ", " control.mode = control.mode, 
##    control.expert = control.expert, ", " control.hazard = control.hazard, 
##    control.lincomb = control.lincomb, ", " control.update = 
##    control.update, control.lp.scale = control.lp.scale, ", " 
##    control.pardiso = control.pardiso, only.hyperparam = only.hyperparam, 
##    ", " inla.call = inla.call, inla.arg = inla.arg, num.threads = 
##    num.threads, ", " keep = keep, working.directory = working.directory, 
##    silent = silent, ", " inla.mode = inla.mode, safe = FALSE, debug = 
##    debug, .parent.frame = .parent.frame)" ) 
## Time used:
##     Pre = 0.407, Running = 9.17, Post = 0.373, Total = 9.95 
## Fixed effects:
##                       mean    sd 0.025quant 0.5quant 0.975quant   mode kld
## (Intercept)         -0.663 0.476     -1.603   -0.661      0.265 -0.661   0
## Fragmentacion        0.171 0.105     -0.035    0.171      0.377  0.171   0
## HogaresNBI          -0.058 0.070     -0.196   -0.058      0.079 -0.058   0
## HogaresHacinamiento  0.120 0.137     -0.148    0.119      0.387  0.119   0
## HogaresJefe.Univ.    0.051 0.022      0.007    0.051      0.095  0.051   0
## Bancos               0.135 0.193     -0.244    0.135      0.513  0.135   0
## 
## Random effects:
##   Name     Model
##     re_u Besags ICAR model
##    re_v IID model
## 
## Model hyperparameters:
##                                                            mean       sd
## zero-probability parameter for zero-inflated poisson_1 5.85e-01 5.00e-02
## Precision for re_u                                     3.20e-01 5.50e-02
## Precision for re_v                                     2.18e+04 1.32e+04
##                                                        0.025quant 0.5quant
## zero-probability parameter for zero-inflated poisson_1      0.481 5.87e-01
## Precision for re_u                                          0.223 3.16e-01
## Precision for re_v                                       7188.279 1.84e+04
##                                                        0.975quant     mode
## zero-probability parameter for zero-inflated poisson_1   6.78e-01 5.94e-01
## Precision for re_u                                       4.41e-01 3.09e-01
## Precision for re_v                                       5.68e+04 1.35e+04
## 
## Marginal log-Likelihood:  -2310.35 
##  is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')

Visualización Espacial

Visualización Interactiva

datos$RR <- res_inflpoi$summary.fitted.values[, "mean"]
datos$RR_LI <- res_inflpoi$summary.fitted.values[, "0.025quant"]
datos$RR_LS <- res_inflpoi$summary.fitted.values[, "0.975quant"]

popup_vars <- c(
  "Fragmentación (Nivel)" = 'Fragmentacion',
  "Valor Tierra Medio ($/m2)" = 'ValorTierra',
  "Hogares NBI (%)" = 'HogaresNBI',
  "Hogares Hacinamiento (%)" = 'HogaresHacinamiento',
  "Hogares Jefe Univ.(%)" = 'HogaresJefe.Univ.',
  "Bancos",
  "Casos",
  "Población" = 'Poblacion',
  "E",
  "SIR",
  "RR_LI",
  "RR",
  "RR_LS"
)
tmap_mode('view')
mapas <-
  tm_basemap(c(Urbano = "OpenStreetMap", Satelite = "Esri.WorldImagery")) +
  tm_shape(datos,
           name = 'Casos')  +
  tm_polygons(
    fill = 'Casos',
    fill.scale = tm_scale_intervals(
      style = "fixed",
      breaks = c(0, 1, 2, 3, 5, 10, 15, 25, 36),
      values = "brewer.yl_or_br"),
    fill.legend = tm_legend(title =  "Casos por Radio",
                            scientific = TRUE,
                            format = "f",
                            digits = 0
    ),
    col = "gray50",
    col_alpha = .5,
    popup.vars = popup_vars
  ) +
  
  tm_shape(datos,
           name = 'Tasa de Infección - SIR')  +
  tm_polygons(
    fill = "SIR",
    fill.scale = tm_scale_intervals(
      style = "fixed",
      breaks = c(0, 5, 10, 15, 20, 25, 30, 50, 60, 90),
      values = "brewer.yl_or_br"),
    fill.legend = tm_legend( 
      title = 'SIR',
      scientific = TRUE,
      format = "f",
      digits = 0
    ),
    col = "gray50",
    col_alpha = .5,
    popup.vars = popup_vars
  ) +
  
  tm_shape(datos,
           name = 'Riesgo Relativo')  +
  tm_polygons(
    fill = "RR",
    fill.scale = tm_scale_intervals(
      style = "fixed",
      breaks = c(0.4, 1, 2, 4, 8, 10, 15, 20, 30, 50, 85, 96),
      values = "brewer.yl_or_br"),
    fill.legend = tm_legend( 
       title = "Riesgo Relativo",
      scientific = TRUE,
      format = "f",
      digits = 1
    ),
    col = "gray50",
    col_alpha = .5,
    popup.vars = popup_vars
  )
mapas
img <-
  "https://ig.conae.unc.edu.ar/wp-content/uploads/sites/68/2022/04/cropped-cromas-68-1.png"

map <-
  tmap_leaflet(mapas) %>%
  leafem::addLogo(img,
                  url = "https://ig.conae.unc.edu.ar/",
                  width = 136,
                  height = 50) %>%
  leaflet::addMiniMap(position = "bottomleft",
                      width = 150,
                      height = 150)
map

# mapshot(map, paste0("Mapa_radios_", format(Sys.time(), "%d_%m_%Y"), ".html"))

Visualización Estática

tmap_mode('plot')

mapacasos <-
  tm_shape(datos)  +
  tm_polygons(
    fill = "Casos",
    fill.scale = tm_scale_intervals(
      style = "fixed",
      breaks = c(0, 1, 2, 3, 5, 10, 15, 25, 36),
      values = "brewer.yl_or_br"),
    fill.legend = tm_legend(title =  "Casos por Radio",
      scientific = TRUE,
      format = "f",
      digits = 0
    ),
    id = "Casos",
    col = "gray50",
    col_alpha = .5) + 
  tm_scalebar(position = c("right", "bottom")) +
  tm_compass(type = "8star", position = c("left", "bottom"))

mapacasos


mapaSIR <-
  tm_shape(datos)  +
  tm_polygons(
    fill = "SIR",
    fill.scale = tm_scale_intervals(
      style = "fixed",
      breaks = c(0, 5, 10, 15, 20, 25, 30, 50, 60, 90),
      values = "brewer.yl_or_br"),
    fill.legend = tm_legend( 
      title = 'SIR',
      scientific = TRUE,
      format = "f",
      digits = 0
    ),
    id = "SIR",
    col = "gray50",
    col_alpha = .5,
  ) + 
  tm_scalebar (position = c("right", "bottom")) +
  tm_compass(type = "8star", position = c("left", "bottom"))

mapaSIR


mapaRR <-
  tm_shape(datos)  +
  tm_polygons(
    fill = "RR",
    fill.scale = tm_scale_intervals(
      style = "fixed",
      breaks = c(0.4, 1, 2, 4, 8, 10, 15, 20, 30, 50, 85, 96),
      values = "brewer.yl_or_br"),
    fill.legend = tm_legend( 
       title = "Riesgo Relativo",
      scientific = TRUE,
      format = "f",
      digits = 1
    ),
    id = "RR",
    col = "gray50",
    col_alpha = .5
  ) + 
  tm_scalebar(position = c("right", "bottom")) +
  tm_compass(type = "8star", position = c("left", "bottom"))

mapaRR

:::